{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vector Spaces: Part I, \"Graphical Representations and Intuition\"\n",
    "\n",
    "\n",
    "Much of statistical learning relies on the ability to represent data observations (measurements) in a space comprising relevant dimensions (features). Often, the number of relevant dimensions is quite small; if you were trying to discern a model that described the area of a rectangle, observations of only two features (length and width) would be all you needed. [Fisher's well-known iris dataset](http://scikit-learn.org/stable/auto_examples/datasets/plot_iris_dataset.html) comprises 150 measurements of only three features 📊 \n",
    "\n",
    "In some cases - particularly with text analysis - the dimensionality of the space can grow much faster. In many approaches to text analysis, the process to get from a text corpus to numerical feature vectors involves a few steps. Just as an exapmle, one way to do this is to: \n",
    "\n",
    "1. **break the corpus into documents** e.g. each on a new line of an input file\n",
    "2. **parse the document into tokens** e.g. split words on whitespace\n",
    "3. **construct a feature vector for each document**\n",
    "\n",
    "One way to accomplish the final step is to consider each token (ie word) as a unique dimension, and the count of each word per document as the magnitude along the corresponding dimension. There are certainly other ways to define each of these steps (and more subtle details to consider within each), but for now, we'll consider this simple one. \n",
    "\n",
    "Using exactly this approach, constructing a vector space from just a few minutes of Tweets (each Tweet considered a document) leads to a space with hundreds of thousands of features! In this high-dimensional vector space, it becomes easy for us to be misled by our intuition for statistical learning approaches in more \"human\" dimensions e.g. one-, two- and three-dimensional spaces. At this point, many people will cite the **\"curse of dimensionality.\"** \n",
    "\n",
    "\n",
    "> There are multiple phenomena referred to by this name in domains such as numerical analysis, sampling, combinatorics, machine learning, data mining, and databases. **The common theme of these problems is that when the dimensionality increases, the volume of the space increases so fast that the available data become sparse.** This sparsity is problematic for any method that requires statistical significance. In order to obtain a statistically sound and reliable result, **the amount of data needed to support the result often grows exponentially with the dimensionality.** Also organizing and searching data often relies on detecting areas where objects form groups with similar properties; **in high dimensional data however all objects appear to be sparse and dissimilar in many ways which prevents common data organization strategies from being efficient.**\n",
    "\n",
    "- [Wikipedia](https://en.wikipedia.org/wiki/Curse_of_dimensionality)\n",
    "\n",
    "This \"curse of dimensionality\" refers to a few related, but distinct challenges with data for statistical learning in high dimensions:\n",
    "\n",
    "- \"increasing dimensions decrease the power of a test statistic\" [ref](http://stats.stackexchange.com/a/30299)\n",
    "- \"Intuition fails us in high dimensions\" [ref](http://homes.cs.washington.edu/~pedrod/papers/cacm12.pdf)\n",
    "- sparse data is increasingly found in the corners/shell of high-dimensional space (this notebook!) \n",
    "\n",
    "I wanted to build more intuition around thinking, visualizing, and generally being more aware of how these phenomena affect our typical kinds of analyses; this notebook is a first step, primarily focused on **building an intuition for inspecting and thinking about ways to inspect spaces when we can longer just plot them.** \n",
    "\n",
    "Along the way, I learned a number of new things, and aim to explore them in follow up pieces.\n",
    "\n",
    "*Note: Beware that there are a lot of reused variable names in this notebook. If you get an unexpected result, or an error, be sure to check that the appropriate data generation step was run!*\n",
    "\n",
    "\n",
    "## Take-aways\n",
    "\n",
    "These are the two high-level objectives we'll aim for:\n",
    "\n",
    "- concepts for inspecting data in high dimensions\n",
    "- illustrations of how high dimensionality squeezes data \"density profile\" to edges\n",
    "    - kmeans `~O(d)`\n",
    "    - \"statistical significance\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import copy\n",
    "try:\n",
    "    import ujson as json\n",
    "except ImportError:\n",
    "    import json    \n",
    "import math\n",
    "import operator\n",
    "import random\n",
    "\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import numpy as np\n",
    "from numpy.linalg import norm as np_norm\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from scipy.spatial import distance as spd\n",
    "import seaborn as sns\n",
    "from sklearn.datasets import make_blobs\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.feature_extraction.text import CountVectorizer\n",
    "from sklearn.feature_extraction.text import TfidfVectorizer\n",
    "\n",
    "sns.set_style('whitegrid')\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple, visualizable spaces\n",
    "\n",
    "We'll start by exploring some approaches to thinking about and inspecting data in spaces that we can comprehend without much effort.\n",
    "\n",
    "## Normal distributions in 2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# number of data points\n",
    "n = 1000\n",
    "\n",
    "# array of continuous values, randomly drawn from standard normal in two dimensions\n",
    "X = np.array(np.random.normal(size=(n,2)))\n",
    "\n",
    "# seaborn plays really nicely with pandas\n",
    "df = pd.DataFrame(X, columns=['x0','x1'])\n",
    "\n",
    "df.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have a 2-dimensional feature space containing 1000 pieces of data. Each coordinate is orthogonal, and we can equivalently think about each data point being represented by a vector from the origin [ (0,0) in 2-dimensional space ], to the point defined by [x0, x1].\n",
    "\n",
    "Since we only have two dimensions, we can look at the bivariate distribution quite easily using `jointplot`. Seaborn also gives us some handy tools for looking at the univariate distributions at the same time 🙌🏼"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sns.jointplot(data=df, x='x0', y='x1', alpha=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another distribution that can provide some hints about the structure of data in a multi-dimensional vector space, is the pairwise inter-point distance distribution for all points in the data. Here's a function that makes this a little cleaner."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def in_sample_dist(X, n):\n",
    "    \"\"\"Create a histogram of pairwise distances in array X, using n bins.\"\"\"\n",
    "    plt.figure(figsize=(15,6))\n",
    "\n",
    "    # use scipy's pairwise distance function for efficiency\n",
    "    plt.hist(spd.pdist(X), bins=n, alpha=0.6)\n",
    "\n",
    "    plt.xlabel('inter-sample distance')\n",
    "    plt.ylabel('count')    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "in_sample_dist(X,n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In unsupervised statistical learning, we're often interested in the existence of \"clusters\" in data. Our intuition in low dimensions can be helpful here. In order to identify and label a grouping of points as being unique from some other grouping of points, there needs to be a similarity or \"sameness\" metric that we can compare. One such measure is simply the distance between all of the points. If a group of points are all qualitatively closer to each other than another group of points, then we might call those two groups unique clusters. \n",
    "\n",
    "If we look at the distribution of inter-point distances above, we see a relatively smooth distribution, suggesting that no group of points is notably closer or further than another other group of points. We'll come back to this idea, shortly. (The inspiration for this approach is found here: [pdf](http://www-users.cs.umn.edu/~kumar/papers/high_dim_clustering_19.pdf))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above, the bivariate pairplot works great for displaying our data when it's in two dimensions, but you can probably imagine that even in just d=3 dimensions, looking at this distribution of data will be really hard. So, I want to create a metric that gives us a feel for where the data is located in the vector space. There are *many* ways to do this. For now, I'm going to consider the **euclidean distance cumulative distribution function\\***. Remember that the [euclidean distance](http://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.euclidean.html#scipy.spatial.distance.euclidean) is the $L_{2}$ norm $dist(p,q) = \\sqrt{ \\sum_{i=1}^{d} (q_{i}-p_{i})^{2} }$ where `d` is the dimensionality of the space. ([Wiki](https://en.wikipedia.org/wiki/Euclidean_distance))\n",
    "\n",
    "\n",
    "\\*in fact, even in the course of developing this notebook, I learned that this is not a terribly great choice. But, hey, you have to start somewhere! `¯\\_(ツ)_/¯` "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def radius(vector):\n",
    "    \"\"\"Calculate the euclidean norm for the given coordinate vector.\"\"\"    \n",
    "    origin = np.zeros(len(vector))\n",
    "    # use scipy's distance functions again! \n",
    "    return spd.euclidean(origin, vector)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# use our function to create a new 'r' column in the dataframe\n",
    "df['r'] = df.apply(radius, axis=1)\n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a couple of ways that I want to visualize this radial distance. **First**, I'd like to see the univariate distribution (from 0 to `max(r)`), and **second**, I'd like to see how much of the data is at a radius less than or equal to a particular value of `r`. To do this, I'll define a plotting function that takes a dataframe as shown above, and returns plots of these two distributions as described.\n",
    "\n",
    "There's a lot of plotting hijinks in this function, so first just look at the output and see if it makes some sense. Then we can come back and dig through the plotting function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def kde_cdf_plot(df, norm=False, vol=False):\n",
    "    \"\"\"Display stacked KDE and CDF plots.\"\"\"\n",
    "    \n",
    "    assert 'r' in df, 'This method only works for dataframes that include a radial distance in an \"r\" column!'\n",
    "    \n",
    "    if norm:\n",
    "        # overwrite df.r with normalized version\n",
    "        df['r'] = df['r'] / max(df['r'])\n",
    "        \n",
    "    fig, (ax1, ax2) = plt.subplots(2,1, \n",
    "                                   sharex=True, \n",
    "                                   figsize=(15,8)\n",
    "                                  )\n",
    "    # subplot 1\n",
    "    sns.distplot(df['r'], \n",
    "                 hist=False, \n",
    "                 rug=True, \n",
    "                 ax=ax1\n",
    "                )\n",
    "    ax1.set_ylabel('KDE')\n",
    "    ax1.set_title('n={} in {}-d space'.format(len(df), df.shape[1] - 1) )\n",
    "\n",
    "    # subplot 2\n",
    "    if vol:\n",
    "        raise NotImplementedError(\"Didn't finish implementing this volume normalization!\")\n",
    "        dim = df.shape[1] - 1\n",
    "        df['r'].apply(lambda x: x**dim).plot(kind='hist', \n",
    "                                               cumulative=True, \n",
    "                                               normed=1, \n",
    "                                               bins=len(df['r']), \n",
    "                                               histtype='step', \n",
    "                                               linewidth=2,\n",
    "                                               ax=ax2\n",
    "                                              )\n",
    "\n",
    "        ax2.set_ylabel('CDF')\n",
    "        plt.xlim(0, .99*max(df['r'])**dim)        \n",
    "        xlab = 'volume fraction'        \n",
    "    else:\n",
    "        df['r'].plot(kind='hist', \n",
    "                       cumulative=True, \n",
    "                       normed=1, \n",
    "                       bins=len(df['r']), \n",
    "                       histtype='step', \n",
    "                       linewidth=2,\n",
    "                       ax=ax2\n",
    "                      )\n",
    "\n",
    "        ax2.set_ylabel('CDF')\n",
    "        plt.xlim(0, .99*max(df['r']))\n",
    "        \n",
    "        xlab = 'radial distance'\n",
    "    if norm:\n",
    "        xlab += ' (%)'\n",
    "    plt.xlabel(xlab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's see these distributions for the 2-dimensional array we created earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "kde_cdf_plot(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a reminder: \n",
    "\n",
    "- the **kernel density estimate** (KDE) is a nice visualization of the \"density profile\", created by assuming there exists a standard normal at each data point, summing all of these curves, and then normalizing the total under-curve area to 1. [The `seaborn` docs](https://stanford.edu/~mwaskom/software/seaborn/tutorial/distributions.html#kernel-density-estimaton) have a nice illustration of this technique. The ticks on the bottom are called a \"rug plot\", and are the values of the data (values of `r`)\n",
    "- the **cumulative distribution function** (CDF) is a measure of the fraction of values which have a value equal to, or lesser than, the specified value, $CDF_{X}(x)=P(X \\le x)$. For the purpose of this session, I want to use this particular to highlight where the observed data is, relative to the \"radius\" of the entire space. The value of the CDF is the fraction of data contained at an equal or lesser \"radius\" value (in d dimensions). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Blobs in 2D\n",
    "\n",
    "Let's add a bit of complexity to the examples above by making the data slightly more irregular: we'll use ``sklearn``'s blob constructor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# data points, dimensions, blob count\n",
    "n = 1000\n",
    "dims = 2\n",
    "blobs = 5\n",
    "\n",
    "# note: default bounding space is +/- 10.0 in each dimension\n",
    "X, y = make_blobs(n_samples=n, n_features=dims, centers=blobs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# convert np arrays to a df, auto-label the columns\n",
    "X_df = pd.DataFrame(X, columns=['x{}'.format(i) for i in range(X.shape[1])])\n",
    "\n",
    "X_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sns.jointplot(data=X_df, x='x0', y='x1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_df['r'] = X_df.apply(radius, axis=1)\n",
    "\n",
    "#X_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time, we'll incorporate one extra kwarg in the `kde_cdf_plot` function: `norm=True` displays the x axis (radial distance) as a fraction of the maximum value. This will helpful when we're comparing spaces of varying radial magnitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "kde_cdf_plot(X_df, norm=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a start, notice that the radius CDF for this data has shifted to the right. At larger `r`, we're closer to the \"edge\" of the space containing our data. The graph will vary with iterations of the data generation, but should consistently be shifted to the right relative to the 0-centered standard normal distribution.\n",
    "\n",
    "Now let's look at the inter-sample distance distribution. Remember that this data is explicitly generated by a mechanism that includes clusters, so we should not see a nice uniform distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "in_sample_dist(X,n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sure enough, we can see that there are in fact some peaks in the inter-sample distance. This makes sense, because we know that the data generation process encoded that exact idea. Since we're intentionally using a data generation process that builds in clusters, we'll always see a peak on the low end of the x axis... each cluster is created with a low (and similar) intra-cluster distance. The other, larger peaks, will illustrate the relationships between the clusters.\n",
    "\n",
    "We may not see precisely the same number of peaks as were specified in the blob creation, though, because we know that sometimes the blobs will be on top of each other and will \"look\" like one cluster. Compare the peaks of this distribution with the pairplot we created with the same data. \n",
    "\n",
    "## Blobs in 3D\n",
    "\n",
    "Let's increase the dimension count by one, to 3, just about the limit of our intuition's abilities. To make the data generation process a bit more reusable, we'll use a function to get the data array and corresponding dataframes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def make_blob_df(n_points=1000, dims=2, blobs=5, bounding_box=(-10.0, 10.0)):\n",
    "    \"\"\"Function to automate the np.array blob => pd.df creation and r calculation.\"\"\"\n",
    "    # nb: default bounding space is +/- 10.0 in each dimension\n",
    "    X, y = make_blobs(n_samples=n_points, n_features=dims, centers=blobs, center_box=bounding_box)\n",
    "\n",
    "    # make a df, auto-label the columns\n",
    "    X_df = pd.DataFrame(X, columns=['x{}'.format(i) for i in range(X.shape[1])])\n",
    "    X_df_no_r = copy.deepcopy(X_df)\n",
    "    \n",
    "    # add a radial distance column\n",
    "    X_df['r'] = X_df.apply(radius, axis=1)\n",
    "\n",
    "    return X, X_df, X_df_no_r, y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n = 1000\n",
    "dims = 3\n",
    "blobs = 5\n",
    "\n",
    "\n",
    "X, X_df, X_df_no_r, y = make_blob_df(n, dims, blobs)\n",
    "\n",
    "X_df.head()\n",
    "#X_df_no_r.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12,7))\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "ax.plot(X_df['x0'],X_df['x1'],X_df['x2'],'o', alpha=0.3)\n",
    "\n",
    "ax.set_xlabel('x0'); ax.set_ylabel('x1'); ax.set_zlabel('x2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sns.pairplot(X_df_no_r, plot_kws=dict(alpha=0.3), diag_kind='kde')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "kde_cdf_plot(X_df, norm=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, compare this CDF to the 2-d case above; note that the data is closer to the \"edge\" of the space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "in_sample_dist(X,n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Higher-dimensional blobs\n",
    "\n",
    "Ok, let's jump out of the space where we can easily visualize the data. Let's now go to d=10. While we can still look at pairwise coordinate locations, we can't see the whole space at once anymore. Now we'll rely on our other plots for intuition of the space profile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n = 1000\n",
    "dims = 10\n",
    "blobs = 5\n",
    "\n",
    "\n",
    "X, X_df, X_df_no_r, y = make_blob_df(n, dims, blobs)\n",
    "\n",
    "X_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# this starts to take a few seconds when d~10\n",
    "sns.pairplot(X_df_no_r, diag_kind='kde', plot_kws=dict(alpha=0.3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "kde_cdf_plot(X_df, norm=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "in_sample_dist(X,n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having seen the way these plots vary individually, let's compare, side-by-side, a similar data generation process (same number of points and clusters) in a range of dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n_points = 1000\n",
    "dim_range = [2, 100, 10000]\n",
    "blob_count = 5\n",
    "\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(2,1, sharex=True, figsize=(15,8))\n",
    "\n",
    "for d in dim_range:\n",
    "    ## data generation    \n",
    "    # random gaussian blobs in d-dims\n",
    "    X, y = make_blobs(n_samples=n_points, n_features=d, centers=blob_count)\n",
    "    ## \n",
    "    \n",
    "    ## calculation\n",
    "    # create a labeled df from X\n",
    "    X_df = pd.DataFrame(X, columns=['x{}'.format(i) for i in range(X.shape[1])])\n",
    "    # add an 'r' column\n",
    "    #X_df_no_r = copy.deepcopy(X_df)\n",
    "    X_df['r'] = X_df.apply(radius, axis=1)\n",
    "    # normalize r value to % of max?\n",
    "    X_df['r'] = X_df['r'] / max(X_df['r'])\n",
    "    ##\n",
    "    \n",
    "    ## plotting\n",
    "    # subplot 1 - KDE\n",
    "    sns.distplot(X_df['r'], \n",
    "                 kde=True,\n",
    "                 hist=False, \n",
    "                 rug=True, \n",
    "                 ax=ax1,\n",
    "                 label='{}-dims'.format(d)\n",
    "                )\n",
    "    \n",
    "    # subplot 2 - CDF\n",
    "    X_df['r'].plot(kind='hist', \n",
    "                   cumulative=True, \n",
    "                   normed=1, \n",
    "                   bins=len(X_df['r']), \n",
    "                   histtype='step', \n",
    "                   linewidth=2,\n",
    "                   ax=ax2\n",
    "                  )\n",
    "    ##\n",
    "    \n",
    "\n",
    "ax1.set_ylabel('KDE')\n",
    "ax1.set_title('n={} in {}-d space'.format(len(X_df), dim_range) )\n",
    "ax2.set_ylabel('CDF')\n",
    "\n",
    "plt.xlim(0, .999*max(X_df['r']))    \n",
    "plt.xlabel('radial distance (%)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, (ax1, ax2, ax3) = plt.subplots(3,1, figsize=(15,9))\n",
    "\n",
    "\n",
    "for i,d in enumerate(dim_range):\n",
    "    X, y = make_blobs(n_samples=n_points, n_features=d, centers=blob_count)\n",
    "    \n",
    "    # loop through the subplots\n",
    "    plt.subplot('31{}'.format(i+1))\n",
    "    # plot the data \n",
    "    plt.hist(spd.pdist(X), bins=n_points, alpha=0.6)\n",
    "    plt.ylabel('count (d={})'.format(d))\n",
    "\n",
    "ax3.set_xlabel('inter-sample distance')  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Text data\n",
    "\n",
    "Most of the time, our unsupervised clustering in high dimensions is a function of using text data as an input. We'll start with a small corpus - again, to build intuition about what the data looks like - and then work up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "small_corpus = [\n",
    "    'The dog likes cats.',\n",
    "    'The blue cat eats brown sharks.',\n",
    "    'Why not, blue?'\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vec = CountVectorizer()\n",
    "\n",
    "X = vec.fit_transform(small_corpus)\n",
    "\n",
    "X.todense()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vec.vocabulary_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's good to remember how to map the matrix-like data onto the words that go into it..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "terms = [x for x,_ in sorted(vec.vocabulary_.items(), key=operator.itemgetter(1))]\n",
    "\n",
    "text_df = pd.DataFrame(X.todense(), columns=terms)\n",
    "\n",
    "text_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "text_df['r'] = text_df.apply(radius, axis=1)\n",
    "\n",
    "text_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "kde_cdf_plot(text_df, norm=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a tiny little corpus, these plots aren't very useful. Let's use a bigger one: this text file (not included in the repo, sorry visitors!) is about a 10-minute, 10% sample of Tweet (body text) from the Firehose. It has a little under 400,000 Tweets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "text_array = []\n",
    "\n",
    "with open('twitter_2016-04-06_2030.jsonl.body.txt', 'r') as infile:\n",
    "    for line in infile:\n",
    "        text_array.append(line.replace('\\n', ' '))\n",
    "        \n",
    "print( len(text_array) )\n",
    "print( text_array[0] )        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vec = CountVectorizer(\n",
    "                    #binary=1,\n",
    "                    ## add dimensionality reduction?\n",
    "                    #stop_words='english',\n",
    "                    #lowercase=True,\n",
    "                    #min_df=10\n",
    "                    )\n",
    "\n",
    "dtm = vec.fit_transform(text_array)\n",
    "\n",
    "dtm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# what fraction of the feature space is full?\n",
    "3051924 / ( 374941*523498 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have to do the radius math slightly differently now, because we're dealing with a scipy CSR matrix instead of a dense numpy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#       (element-wise sq) (row sum) (flatten) (sqrt)\n",
    "dtm_r = dtm.multiply(dtm).sum(axis=1).A1**0.5\n",
    "\n",
    "\n",
    "#print(len(dtm_r))\n",
    "#print(dtm_r)\n",
    "#print(min(dtm_r), np.median(dtm_r), max(dtm_r))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "s = pd.Series(dtm_r)\n",
    "\n",
    "plt.figure(figsize=(15,6))\n",
    "s.plot(kind='hist', \n",
    "       cumulative=True, \n",
    "       normed=1, \n",
    "       bins=len(dtm_r), \n",
    "       histtype='step', \n",
    "       linewidth=2\n",
    "      )\n",
    "\n",
    "plt.ylabel('CDF')\n",
    "#plt.xlim(0, .99*max(dtm_r))\n",
    "plt.xlim(0, 6)\n",
    "plt.xlabel('radial distance')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# This is a super interesting side note: some tweets can totally throw off your distribution. \n",
    "# This one Tweet had 114 repetitions of a single character. If you swap the xlim() commands \n",
    "#  above, you'll see that r extends to over 100. This is why:\n",
    "#text_array[ s[ s > 114 ].index[0] ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `<record-stopping screeching noise>` \n",
    "\n",
    "Ok, so I spent some time working with this data, and I'll be honest: I expected this distribution to be much more skewed to large `r`! In fact, I thought it would be *more* exaggerated than the blob examples above. \n",
    "\n",
    "Since I didn't have enough time to dig any deeper for this session, let's keep this observation in the back of our minds, and come back to it in another session. \n",
    "\n",
    "We can round out today's discussion with one more relevant topic..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dimensionality reduction\n",
    "\n",
    "Before we end this session, we'll consider one more facet of high-dimensional spaces: reducing them to lower dimension. For now, we'll illustrate the effect of using principal component analysis using the same inspection techniques we've been using all along.\n",
    "\n",
    "If we try to densify the 500k+ dimension document term matrix above, we'll run out of RAM. So, let's use a synthetic data set. \n",
    "\n",
    "First, we look at our metrics in 10,000 dimensions, then after PCA to bring them down to 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n = 2000\n",
    "dims = 10000\n",
    "blobs = 10\n",
    "\n",
    "\n",
    "X, X_df, X_df_no_r, y = make_blob_df(n, dims, blobs)\n",
    "\n",
    "#X_df_no_r.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "kde_cdf_plot(X_df, norm=True)\n",
    "plt.xlim(0,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "in_sample_dist(X,n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we know that the data generation process built in the notion of identifiable clusters. Let's see if we can surface that information by projecting our high-dimensional data and space down into a smaller number using principal component analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# now apply PCA and reduce the dimension down to 3\n",
    "pca = PCA(n_components=3)\n",
    "\n",
    "X_df_3d = pd.DataFrame(pca.fit_transform(X_df_no_r), columns=['x0','x1','x2'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# add in that radial distance column\n",
    "X_df_3d['r'] = X_df_3d.apply(radius, axis=1)\n",
    "\n",
    "X_df_3d.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# add in the labels so we can color by them\n",
    "X_df_3d['y'] = y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# nb: using the vars kwarg seems to remove the ability to include KDE\n",
    "sns.pairplot(X_df_3d, \n",
    "             vars=['x0','x1','x2'], \n",
    "             hue='y', \n",
    "             palette=\"husl\",\n",
    "             diag_kind='kde',  \n",
    "             plot_kws=dict(s=50, alpha=0.7)\n",
    "            )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "kde_cdf_plot(X_df_3d, norm=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#in_sample_dist(X_df_3d[['x0','x1','x2']],n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given the two plots just above, it seems like we've both done a good job of representing the underlying clusters in our lower-dimensional space, and moved the data away from the extreme edges of the feature space. We should expect that both our algorithms can run more efficiently (faster), and achieve a higher level of significance.\n",
    "\n",
    "## Still to come...\n",
    "\n",
    "In future installments, I look forward to:\n",
    "\n",
    "- a deeper understanding of our typical high-dimensional text feature space \n",
    "- strategies for dealing with distance calculations in high dimensions\n",
    "- ... probably some other stuff, too..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
